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Abstract 

We consider a continuum model describing the dynamic behavior of nematic 
liquid crystal elastomers (LCEs) and implement a numerical scheme to solve 
the governing equations. In the model, the Helmholtz free energy and Rayleigh 
dissipation are used, within a Lagrangian framework, to obtain the equations 
of motion. The free energy consists of both elastic and liquid crystalline con- 
tributions, each of which is a function of the material displacement and the 
orientational order parameter. The model gives dynamics for the material dis- 
placement, the scalar order parameter and the nematic director, the latter two 
of which correspond to the orientational order parameter tensor. Our simula- 
tions are carried out by solving the governing equations using an implicit- explicit 
scheme and the Chebyshev polynomial method. The simulations show that the 
model can successfully capture the shape changing dynamics of LCEs that have 
been observed in experiments, and also track the evolution of the order parameter 
tensor. 



1 Introduction 

Liquid crystal elastomers (LCEs) are orientationally ordered solids, combining 
features of liquid crystals and elastic solids. They were first proposed by de 
Gennes [5] and first synthesized by Finkelmann et al [5] . They consist of weakly 
cross-linked liquid crystal polymers with orientationally ordered side or main- 
chain mesogenic units. They exhibit many new phenomena not found in either 
liquid crystals or polymers. The salient feature of LCEs is the strong coupling 
between mechanical deformation and orientational order. As a consequence of 
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this coupling, mechanical strains change the order parameter and hence physical 
properties of LCEs, and, conversely, extermal stimuli, such as light, affecting 
orientational order can produce large shape changes [H [T31 [13 HH] • 

Although many fascinating experimental results have been obtained studying 
the dynamic response of LCEs to external stimuh PlHll^fn i fTi l fTT l l^l^. 
their dynamics is not fully understood. In this paper, we implement a non-local 
continuum model [3, chose and explicitly define a specific representation and 
carry out numerical simulations to explore the dynamic behavior. Our work 
thus includes both components of modeling and simulation. 

In the model, the Helmholtz free energy and Rayleigh dissipation are com- 
bined, using a Lagrangian approach, to obtain the dynamics. The free energy 
consists of both elastic and nematic contributions, and includes volume conser- 
vation. As a special case of the continuum model, we choose a simple local form 
of the nematic free energy, the Maier-Sauper free energy, to describe nematic 
contributions. Our model considers only the uniaxial phase of nematic LCEs, 
for which the order parameter tensor can be expressed in terms of a scalar order 
parameter and nematic director; direct contributions to the free energy from 
spatial variations of the order parameter and director are neglected. These sim- 
plifications make our model more tractable both theoretically and numerically. 
Subsequently, the governing equations can be derived explicitly using both con- 
served and non-conserved order parameter dynamics. We thus obtain the time 
dependent equations for the displacement, scalar order parameter and nematic 
director. 

The equations obtained are more complicated than the standard Navier- 
Stokes equations in the Eulerian frame. First, besides the pressure term and 
the viscous term, there is also an elastic term in the velocity equation. Second, 
our equations are written in a Lagrangian frame. This choice is straightforward 
for capturing the orbit as well as the dynamics of each particle in the LCE 
sample. Indeed, Eulerian coordinates are not well suited to our problem since 
the domain occupied by the LCE sample varies in time. Moreover, the derived 
velocity equation is very stiff due to the presence of different time scales in 
the problem, posing challenges for the simulation. The simulation is therefore 
a fascinating but a very formidable problem. In this work, we employ the 
Chebyshev polynomial method 15^ to discretize the spatial derivatives in the 
dynamical equations. This method, as a typical spectral method, can achieve 
high accuracy, and is particularly well suited to our simulation as our system is 
non-periodic. We also apply the popular implicit-explicit (IMEX) schemes for 
the time-discretization of the equations. Specifically, a combination of second- 
order Adams-Bashforth method for explicit terms and Crank-Nicolson method 
for implicit terms [3 [TUl [H] are used. 

The paper is organized as follows. In Section 2, we give the details of the 
model as well as the derivation of the governing equations. To obtain the equa- 
tions, we first calculate the functional derives of those functionals with respect to 
the principal variables, i.e., material displacement, order parameter and nematic 
director, and then apply the appropriate conserved/non-conserved dynamics for 
each of the variables. In Section 3, we present the numerics for solving the 
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equations and then show the results of our simulations. Conclusions are given 
in Section 4. 



2 Modeling nematic LCEs 



To describe the dynamics of LCEs, in addition to orientational order, one needs 
to track the time evolution of the position of the crosslinks of the LCE network. 
In the case of uniaxial nematic LCEs, the sample can be characterized by the 
displacement, order parameter and nematic direction at each Lagrangian lattice 
site, corresponding to a cross-link. Our work is to study how these key variables 
evolve in time when the sample is subjected to external stimuli. To this end, we 
represent the continuum model in terms of these variables, derive the governing 
equations and implement the simulation by solving the equations numerically. 

Our continuum model consists of the elastic free energy density and the ne- 
matic free energy density with coupling between orientational order and defor- 
mation of the network, a Rayleigh dissipation function and a volume preserving 
functional. The governing equations are derived from these by applying the 
appropriate dynamics for each key variable. In what follows, the functionals 
and the functional derivatives are discussed. 

2.1 Energy functionals 

The free energy in our model is composed of elastic and nematic contribu- 
tions. The elastic free energy describes the nonlocal interaction between con- 
nected cross-links of elastomers, while the nematic free energy represents the 
anisotropic dispersion interactions of the mesogenic constituents. 

2.1.1 The elastic free energy 

Let ot denote a material point in the LCEs sample, and x(Q:,t) the location of 
the point a at time t. To describe the nematic ordering in the LCEs crossing- 
link network, an effective dimensionless step length tensor L is introduced, being 
written as L = I -I- where Q is the orientational order parameter tensor, 

I is the identity matrix and fi is the relaxation parameter. In the uniaxial 
phase case, Q = ^dnn^ — ^1), where n represents the unit vector along the 
average alignment direction of the molecular symmetry axes, S is the scalar 
order parameter describing the degree of alignment of the molecular axes with 



In the undeformed state, the probability density of finding in the LCEs 
sample a polymer chain of length C starting at a and ending at a.' can be 
written as 
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where Lq = I + 2^Qo is the effective step length at the initial state [15] . Since 
Qo is assumed to be slowly varying compared to the distance between crosslinks, 
we evaluate Qo and Lq at position a. 

At time t, the probability density of finding a polymer ending at x(a,t) 
and x(a',t) shares the same form as Po{a.,a') with Lo(a) being replaced by 
L(Q:,t). The free energy of the particular polymer initially ending at a and 
a' is —kT\ii[P{x{a.,t),x{a',t))], where k is the Boltzmann's constant, T is the 
temperature. The total elastic free energy at time t is 



el 



= ^J d^a J d^a' pcPo{a,a')(~kT\n[P{x{a,t),x{a' ,t))f^ 
= Jd^aJ d3a'i?(a,a')(^Wa',i)-x(a,t))^L-i(x(a',t)-x(a,t)) 
"^IndetL), (1) 
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where J-ei is the elastic free energy density, 

and pc is the number density of crosslinks. 
2.1.2 Nematic free energy 

Perhaps the most successful description of nematic order is Maier-Saupe theory. 
Here, the single particle potential is 

£ = -UpicSP2{cose) + i[/p/c52, 

where U is an interaction strength, pic the number density of the liquid crys- 
talline constituent, and S the scalar parameter defined as S =< P2{cos9) >, 
where 9 is the angle between the symmetry axis of a mesogen and the nematic 
director n. 

The nematic free energy can be written as: 



Fne,n = J O. (^pickT\B.{ j exp{--^)dn) 



I d a. {■:^PiJJS - pickT\n{ / exp( — )dn 



(2) 



where dfi = sm6d9d(j>, and 9 is the polar while (j) is the azimuthal angle. 

To simplify the expression for Fnemi we take the Taylor's series expansion of 
the integrand, and obtain a Landau-de Gennes form for the free energy density: 

-Fnem = \CaS^ ~ ^^5^ + ^C^S" + 0(5=^), (3) 
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where 

In the expressions for Ca and C^, hkTjpi^ = T /T* h/Q , and p;c = 
so we write Ca and Cb as 

= 500(|J- - l)pcfcT, 

Cb = AQQpckT. 
Ca = bOOpckT 

where T* ~ 355K, the hmit of undcrcooUng of the isotropic phase, is very near 
the nemati-isotropic transition temperature. 

The nematic free energy can therefore be approximated by 



2.1.3 The Rayleigh dissipation function 

The total dissipated power in the system is 



J d^ocTZa (4) 



and the Rayleigh dissipation function (dissipated power /volume) is 

7^d = ^7iD : D + 72D : Q + ^73(3 : Q, 

where D = (VxU + VxU^)/2 is the symmetric rate-of-strain tensor, u = x and 
7i, i = 1, 2, 3 are viscosities. 

To be consistent with the variable of integration cx in the expressions for 
the free energy, we rewrite the rate-of-strain D in terms of the Lagrangian 
coordinates. Indeed, note that the relation 

Vax = (Vxu)F, 

where F = ^ is the deformation gradient, then 

D = i [(Vax)F-i +F-^(Vax^)] • 
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2.1.4 Volume preserving functional 

The above free energy presents no restrictions on tlie sample volume. It is 
known, however, from experiments, that most rubbers and LCEs are nearly 
volume conserving [TTl [IB] . We therefore introduce a term controlling volume: 

where J = det(F) and A is a positive constant. 



2.2 Derivation of the governing equations 

The equations of motion are determined via a Lagrangian approach, by ex- 
tremizing the action in the presence of dissipation. The Lagrangian is given 

by 

where Skin is the kinetic energy density, and J- is the free energy density. The 
equations of motion are given by 

/ d^'ai-— - — + — = (5 
J at ox ox ox 

/ d^a{- r - tt; + ^) = (6) 

J ^dtdS SS 5S' 

and 

o d (9£ 6C 6TZ. ^ 
dt on on on 



To derive the equations of motion for the displacement (x) , order parameter 
(5) and nematic director (n), one needs to calculate the functional directives 
of Fei^ Fnerrn F^oi with rcspcct to X, S, n, and of the kinetic energy Ekin and 
dissipation TZ. with respect to x, S and n. 



2.2.1 Functional derivatives of elastic free energy 

The elastic free energy is given by 

Fei{a) = j d^OLFel{a) 

= j d^ad^a'H{a,a')(^-^{x{a',t)-x{a,t))'^L-^{x{a',t)~x{a,t)) 
+ ^lndetL) (8) 

We calculate the functional derivatives by the standard procedure, which we 
illustrate here in computing SFei/Sx. Let y — y{cy.) be an arbitrary function. 
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For convenience, we denote = Xi(Q:',<) — Xi{a.,t), hi = yi(a') — yi{a.), and 
g{'x.) = (x(a', t) — x(a, t))'^'L{a., t)~ (x(q:', t) — x(a, t)). Then for any e > 0, 

5(x + ey) = Lr.i(a,t)(a, + eb,)(a, +ebj) 



One gets 



lirrie 



g(x) + eL-.i(a,t)(a,b, + a.b,) + 0{e^). 



,o^(5(x + ey) = L^/(Q;,i)(a»bj +ajb,), 



and then 



i^ei(x + ey)|e=o = y dPa j d^a' H {a, a') {^)Li/ {a,t){a,h ^ + a^h,) 



-L,/(a',t)]aiyj(a) 



Therefore, the functional derivative of the elastic free energy F^i is 



SF, 



X J 2Cb 



(9) 



We proceed similarly to calculate SFei/SS, and get 



SS 



= J d^a'H{a,a') 



1 



2(l-/xS')(l + 2^S') 



2/:fe (1 - ^^5)2 



(x(a) - x(a )) I 



(2^5' 



±^nn^) (x(a)-x(«')) 

(10) ■ 



SF, 
S: 



We next find SFei/Sn. Proceeding as above, we obtain 

—6fiS 



(n- (x(a) -x(a')))^n 



(n.(x(a)-x(a')))(x(a)-x(a')) 

(11) 



The purpose of the non-local description in the model rather than a gradient 
expansion is to ensure that no artifacts arise in the dynamical equations due 
to truncation of the gradient expansion when the variation is carried out [T2] . 
Once the functional derivatives have been evaluated, gradient expansions can 
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safely be carried out. Since evaluation of the integrals is cumbersome and 
computationally expensive, we therefore turn here long-wavelength expansions 
of the integrands, and and apply these to SF^i / 6x, SF^i /SS and 6Fei / 6n. 
Let a' - a = {^y^/^hl^^fS, then 



and 



rf-''«'^(^)-V2(det(Lo)i/2rf-^/3, 
H{cx, a') = (ipfcT)( 2^)3/2(det Lq)" exp(-/3 • /3), 



To make the expression of dF^i/Sx clear, we here consider its i*'' component 



i^-^)i = -'^^-'^'ilpckT) I d^/3exp(-/3./3)[Ly-i(a)^|j^Lj/|„/3„Lif>„ 



(12) 



We recall the useful identities: if Mp = rf exp{—r]^)dr], then Mq = y/n, 
M2 = Mo/2, and M4 = 3Mo/4. Then °° 



or 



^ = -(ipefcT)Va-(L-^FLo). 



(13) 
(14) 



Proceeding similarly, 

SF,i 
SS 



7r-3/2(ip,fcr) 1 d3/3exp(-^./3)^^— 



3(l + 2M^g^) ^ 1/2^ 
" (1 + 2^5)2 ^° ^ 

+ ^-'/'(^PcfcT) y d3/3exp(-/3 



/3) 



-3^25* 



(l-/x5)(l + 2/i5) 



~ ^2^^^^^2(1-m5)2 



tr(FLoF^) - fc|^tr(nn^FLoF^) 



(1 + 2/i5)2 



_ 1 Sti'^S 



(15) 
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and 



S: 



p , 1 r 



— 6/i5 



m5)(2m5+1) 



(n^FLy^/3)2n 



(n^FLy^/3)FLy^/3 
(16) 



For simplicity, we denote G = FLy^ The pth element of / d^/3exp{-\/3\^){n'^Gl3)G/3 



IS 



^njGjjGpfc(5jfc = ^GpfeG^.nj. Therefore 



/3/2 3/2 
d3/3exp(-/3 • /3)(n^G/3)G/3 = ^GG^n = !L_FLoF^n. 

Moreover, as (n'^FLy^/3)2 = /S^L^^F'^nn'^FLy^/S, one gets 

3/2 



/ o / 2 
d^l3exp{-l3 ■ /3)(n^FLy^/3)2 = ^tr(nn^FLoF 



Consequently, 



^(l-^5)(l + 2/i5) 



[FLoF"^ - tr(nn^FLoJ^'^)I]n. (17) 



2.2.2 Functional derivatives of nematic free energy 

As the nematic free energy can be approximated by 

The functional derivative of the free energy density with respect to S is 



ss 



CaS — CffS + CcS , 



and the functional derivatives with respect to x and n are all zero. 

2.2.3 Functional derivatives of volume preservation functional 

The volume conserving term does not depend on x, S, n, S and h. We want 
to find SFyoi/Sx. 

Proceeding as before, we obtain 



j;o;(x + ey) 



o=Ay d3a[tr(F-i|^)(J-l)J 



(18) 



and, after integrating by parts and requiring J = 1 on the boundary, we obtain 



Jx 



= -Va-(A(J-l)JF-^). 



(19) 



2.2.4 Functional derivatives of Rayleigh dissipation 

We write the dissipation as 



7^ 



1 



d-^a -7iD : D + 72D : Q + -73Q : Q = 7^l + 7^2 + 7^3 (20) 



Since TZi is independent of x, S, n, S and ri, we calculate the functional 
derivative 6TZi/6x. Proceeding as before, we obtain 



d_ 



=0 = 11 / d^a 



tr(DF-^Vay^) 



and integration by parts gives 



-71 Va • (DF-^). 



(21) 



A calculation similar to that above shows that 

6112 



-72 Va • (QF-^ ). 



We note that Q = ^dnn^ — |l) + |6'(hn^+nh^). Then simple calculations 
yield 

= 72D:(-nn - -I 



and 



^5 
6h 



= 3725 



Dn - (n^Dn)n 



Similarly, we obtain 



and 



'57^3 3 . 



<57^3 9 2 . 
7h = 



By combining these, we obtain the functional derivatives of the Rayleigh 
dissipation as 



^ = -71 Va • (DF-^) - 72 Va • (QF"^), 
ox 



(22) 
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5n ^ f'i T U 

— :- = 72D : -nn — -I 

SS \2 2 



and 



^ = 3725[Dn- (n^Dn)n] + 



(23) 



(24) 



2.2.5 The equations of motion 

We now derive the equations of motion governing the time evolution of the 
velocity u of the elastomer, and of the nematic order parameter Q, expressed 
in terms of S and n. 
Recall that 

jC=j d^Ci{£un - J") 
where Skin is the kinetic energy density, T is the free energy density, and 



d-'oL — — - — —\ -0 

dt &k Sx S-k J 

d dC 5C Sn 

dt ds 5s 6s 



and 



d dc _sc sn 

dt dh 6n 6h 



= 0. 



(25) 

(26) 

(27) 



First the Lagrangian map x{a.,t). Recasting the kinetic energy from the 
current to the initial configuration gives: 

^kin ~ I d xSf^in 

Jn{t) 



= [ d^x[l-pm{x,t)u{x,t) ■U{x,t) 
Jn{t) \ ^ y 

= / rf'^a ( ^p„(x(a,t),i)x(Q;,t) •x(a,i)J(a,t) ] 

JS2(0) / 

= \ d^ct \-pmici, 0)x(a, t) ■ x(a, t)\ , (28) 
Jn(o) / 

where we use mass conservation pm(x(a, i), t) J(a, i) = ^^(0,0), with the 
mass density of the hquid crystal elastomer. We then have 

= Pm^ = PmU, (29) 

and subsequently 

d . , 6{F^i + Kern + F^oi) ^STZ ^ 

— (p„u) + + — = 0, 

at ox ox 
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or 

Pr,^^ = {^PckT)Va ■ (L-^FLo) + Va ■ (A(J - 1) JF"^) 

+ 7iVa • (DF-^) + 72Va • (QF-^). (30) 

We remark that the couphng of strain and orientational order, the sahent aspect 
of hquid crystal elastomers, is implicit in the first term of the RHS. Since L = 
I + 2/iQ, spatial variations of the order parameter give rise to stresses, and in 
turn, to elastomer motion. 

Next, consider the dynamics of the order parameter expressed through the 
variables S and n. Since the Lagrangian C does not depend on S or ri, the 
equations of motion give 

d{Fel + Fnern + F^ol) , dTZ _ ^ 



dS dS 



dn dh 

or 



0, 



2 dt '2^' '2{1- iiS)^ 



^ 72D : gnn^ - il) , (31) 

and 

Eqs. ([50]) . (PT|) . and are the equations of motion for the nematic liquid 
crystal elastomer system. 

We make these equations of motion nondimensional by introducing the fol- 
lowing dimensionless quantities: 

, u , q; ft 

u = -, = T' ^ = 

U T 

where b is the step length of liquid crystal, and constants u and t are to be 
determined. 

The equation ([50)) becomes 

P™^^ = ilpckT)^^a' ■ (L-^FLo) + ■ (A(J - 1) JF-^)+ 

+ 7i^Va. • ((Vc,'uF-i + F-^Va'U^)F-^) , (33) 
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where we let 72 = for simplicity for the time being. Letting u = b/r and 
T = Isl PckT^ the above equation reads: 



= \^^' ■ (L^'FLo) + Va, • (A'( J - 1) JF"^) 
+ ^^OL' ■ ((V^'uF-i + F-^Va'U^)F-^ 

where A = PchTpmii^/^l and A' = K/pckT. 
With this choice of parameters we have 



(34) 



1 



dS_ 



tr(FLoF^) - 4^±^4^tr(nn^FLoF^) 



(1 + 2fiSy 



200 / 

{1 - fiS){l + 2iIS) ^ ~3~V 



-5( 



Tern 
"360" 



1)S + AS^~5S'' 



(35) 



where Tern — Tem(Q;',t') is a temperature function depending on location a' 
and time t' . 



pS 



dn 



dt' 352(1 -^5)(1 + 2/^5*) 



FLoF^ - tr(nn^FLoF^)I 



(36) 



This equation preserves the length of the director n, as required. 
Finally, the deformation matrix F satisfies 

while (again) the Lagrangian map x satisfies 

9f 



(37) 



(38) 



3 Numerical results 

We present simulations of the dynamics of an LCE sample - using Equations 
([34l - [38|l ~ when exposed to external illumination and subject to two different 
boundary conditions. The LCE sample is taken as box-shaped, as in Fig. [T^. 
In the first set of simulations, zero-stress boundary conditions are imposed over 
the sample surface (i.e., the sample is "free"). In the second case, one end 
of the sample is anchored to a wall, with the remainder free. In either case, 
gravitational loads are neglected. Numerically, the difference between these 
two cases lies only in the treatment of the velocity on one face of the sample. 
However, the dynamics the two cases present is quite different, as observed in 
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3.1 Methods 



To discretize the equations of motion, one needs to consider suitable schemes 
for approximating both spatial derivatives and time derivatives, here we em- 
ploy the spectral Chebyshev polynomial method to discretize spatial deriva- 
tives with high efhciency and accuracy [TS]. As for the time-discretization, we 
use a popular implicit-explicit scheme that is a combination of second-order 
Adams-Bashforth scheme for the explicit term and Crank-Nicolson scheme for 
the implicit term P HTUIITS] . 

We now outline the Chebyshev polynomial method and the implicit-explicit 
time-stepping method. The dynamics is simulated in the Lagrangian domain 
ri(0) — [—a, a] X [— fo, 6] x [— c, c], which by definition is fixed in time. This is 
trivially mapped to the cube [—1,1]^. This cubic domain is then discretized 
in each direction on the Gauss-Lobatto points (e.g. in the first coordinate, 
cxij = cos{jTr/N), j — O(l)A^). This allows allows spatially dependent fields, 
such as u or F, that are represented discretely on these points to also be repre- 
sented efhciently, via FFT, as finite sums of Chebychev polynomials [TS]. The 
Chebychev representation can then be used to provide highly accurate deriva- 
tive approximations upon the grid. To illustrate in one dimension, let u{x) be 
defined on [—1, 1] and approximated by un{x) = X^^o ^pTp{x) where Tp is the 
p*'' Chebychev polynomial. The Cp's are determined by requiring M^r to inter- 
polate u at the Gauss-Lobatto points. This Chebychev representation allows us 
to construct approximations to u'^p\x), at the Gauss-Lobatto points, that have 
the form 

N 

UN\x^) = Ww(a^J-), « = 0(l)iV, 

j=0 

where [d^^^] is the the p*'' Chebyshev differentiation matrix (see [H]). Since our 
three-dimensional grid is of tensor product form, derivatives are easily gotten by 
application of such matrices along lines of constant coordinate of the discretized 
data. 

The LCE dynamics in which we are interested takes place in the over-damped 
regime, that is the "Reynolds number", A73/71, associated with viscous fluid 
damping is very small. Hence, if we retain Ut in the dynamics, we must implicitly 
treat the viscous damping so as to avoid extreme constraints on the time-step 
that would be imposed by using an explicit scheme. Here, we choose a pop- 
ular implicit-explicit method which is described for a typical time dependent 
equation: 

^ = /(w) + i^5(u), 

with /(•), g{-) being nonlinear and linear terms, respectively. We apply a second- 
order Adams-Bashforth method to the nonlinear terms, and Crank-Nicholson 
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averaging to the linear term, or 

n+l _ ,,n o 1 ,. 

^f— = 2/K) - ^fiu"-') + 2 + (39) 

where At is the time step size and m" is the approximation of u{n At). This 
scheme involves solution values on three time levels. The first time-step is taken 
by setting = = u{0) (see [15]'). 



3.2 Treatment of interior and boundary points 



Note that the equations ((351) ^^nd (|36)) involve no spatial derivatives, and so are 
solved directly by the second-order Adams-Bashforth method and without using 
the Chebyshev approximation. More care must be taken with the momentum 
equation (|34p as it is a source of stiffness in the numerical treatment, and since 
its advancement involves boundary conditions. This is an important issue, since 
distinct boundary conditions result in completely different behaviors of the LCE 
sample. In what follows, we mainly discuss how to handle the equation for both 
boundary and interior points. 

The three-dimensional grid using the Gauss-Labatto points is composed of 
both surface and interior points. To update the velocity at the interior points, 
we need to solve a large linear system gotten by applying the implicit-explicit 
method for the velocity equation ([34]), coupling this to boundary conditions. 
For a "free" LCE sample, this is a condition of zero normal stress, which under 
discretization of velocity gradients in the viscous stress provides a coupling 
condition of the boundary velocities to the interior velocities. 

The velocity equation (IMl) can be discretized as follows: 
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u" 



At 



3 
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1 

2 
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473 
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1 



(A(J-l)JF-^ 
{A{J-l)JF-^- 
(Vu"+i(F-i)"H 
(Vu"(F-i)"-i H 



^L-iFLo) 



(F-^)"(Vu^)"+i(F-^)" 
(F-'r)"^i(Vu^)"(F-^)" 



where [•]" represents the value at nAt, for instance, (F~-^)" denotes the value of 
F~"^ at nAt. This equation amounts to a large linear system for the unknown 
velocity u""*"^ at the interior points. Couplings within the matrix arise through 
expansion of spatial derivatives (that is, gradients and a tensor divergence) 
through the Chebyshev expansion of the velocity. Despite the many entries in 
the matrix generated through the derivatives, the matrix is nonetheless rather 
sparse, and we explicitly construct the entries and store this sparse matrix. 
Once the boundary conditions are appropriately integrated, we solve this large 
system using the iterative GMRES method [16] . 

Surface values of velocity are either additional unknowns, or are specified as 
in the case of having an anchored surface where u = on that face. The former 
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is the case of the zero stress boundary condition. In the Lagrangian frame, using 
Nanson's formula [T9l, this boundary condition can be written as 

A(J - 1)1 + iL-iFLoF^l • JF-^i^o + ^ [(VauP-i 
L 2 J 273 L 

+F-^VaU^)] • JF-^iya = 0, (40) 

where vq denotes the outward normal unit vector to the time invariant sur- 
face 9J7(0) = dil. In our method, this vector 1^0 takes different values at 
the side, edge and corner points on the boundary of the cubic volume of the 
LCE sample. Specifically, at the side points, 1^0 takes value from the set 
{(±1, 0, 0), (0, ±1, 0), (0, 0, ±1)}; at the edge points, the set becomes {^(±1, ±1, 0), 
^(±l,Tl,0),^(±l,0,±l),^(±l,0,Tl),^(0,±l,±l),^(0,±l,Tl)};and 
at the corner points, the set is then {-^(±1,±1,±1), -^(±1,±1,=f1), -^(=f1,±1,±1), 

;^(Tl,±l,Tl)}. 

We rewrite the boundary condition equation (|40p as: 



[(VauF-i + F-^VaU^)] ¥'^1^0 = Ai^o, 

where A = -^[iL^^FLo + A(J - 1)JF"^]. This boundary condition is 
evaluated at the {n + l)st time-step. Given that all of the dynamics equations, 
bar that for the velocity u, are treated explicitly, we can consider F and A as 
being considered known, and the velocity gradients as unknowns. Gradients are 
either tangential to the sample surface, and hence only couple together boundary 
points (upon approximation of gradients using the Ghebyshev representation) , 
or normal to the surface and hence couple together boundary and interior points. 

This yields a closed set of equations for u"+^ at the interior points, and on 
those surface upon which a zero stress boundary condition is imposed. As said 
above, this system is solved via the GMRES iterative method [T5] . 

An important issue that should be emphasized is how to employ the GMRES 
method efficiently. Indeed, note that the coefficient matrix is very sparse, we 
just need to store the non-zero entries of the coefficient matrix, and also the 
row and column indices of these non-zero entries. Then in the matrix- vector 
multiplication, these non-zero entries will be multiplied by those elements of the 
vector using the above stored column indices. This procedure can save lots of 
memory and also accelerate the matrix- vector multiplication considerably. 



3.3 Simulations 

Before proceeding, we discuss the choice of dimesionless parameters. These in- 
clude the coefhcient of acceleration A — PchTp^b'^ /'^i, the viscosity ratio 71/73, 
the coefficient for volume conservation A, and the anisotropy of step length /i 
appearing in the tensor L. Taking typical values [13] we have A = O{10~^) 
and 71/73 — O(10^~^). Hence, inertial forces in the material are quite small. 
Ideally, we should choose a very large value of A to enforce material incompress- 
ibility, but this imposes a severe time-step restriction in our numerical scheme; 
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We use A ~ 10^. The parameter lies in the range [0, 1.0]. A large value of 
H corresponds to a large order parameter, which accelerates the deformation 
process of the LCE sample. We use /j, = 0.9 in the simulations. 

We now consider the simulated dynamics of the first case of a "free" LCE 
sample being exposed to illumination from above. In this simulation, the sample 
size is 8 X 8 X 1, with Ni = 32, N2 = 32 and N3 = 10 points being used in each 
direction, respectively. The initial data used was Uq = 0, ng = y, Xq = a, and 
So = s, where s is the constant value found as the minimizer of the Landau- 
de Gennes free energy density (O given a uniform temperature throughout the 
sample corresponding to 290K (in dimensional units). As discussed earlier, we 
neglect thermal diffusion and assume that the temperature is uniform in each 
horizontal slice of the sample, decreasing linearly from top (420K) to bottom 
(290K). 

Figure [T] shows the deformation process from this initial configuration. The 
final result - a saddle shape - is very similar to that observed in the experi- 
ment of Palffy-Muhoray et al. (see Fig. 4 of [13] )■ The evolution proceeds in 
three stages: an initially slow and small bending, followed by rapid and large 
deformation, and finally a slow relaxation to a terminal shape. This dynamics 
is driven by the evolution of the orientational order parameter, s, as it adjusts 
its values (low on the top and higher on the bottom) in response to the imposed 
temperature gradient. The inhomogeneous spatial distribution of orientational 
order, especially through the thickness of the LCE sample, gives rise to large 
stresses and hence creates a strong driving force towards changing the shape of 
the sample. 

The last two plots of Fig. [1] show the late-time deformed sample from two 
difl[erent perspectives. Here one finds that the length along the y— axis has 
become shorter, while that along the cc— axis has increased. This is again due to 
the time evolution of order parameter. At the top surface, given its increased 
temperature, the degree of order of the rod-like mesogens decreases. Since 
these mesogens are initially aligned along the y— direction, this loss of order 
leads to a contraction of the sample along the y— direction and corresponding 
extensions along the x— and z— directions. Since the temperature is different on 
each horizontal layer, the degree of contraction also is also different. It is this 
difference in contraction and expansion through the thickness of LCE sample 
that results in the observed saddle-shaped deformation. 

The simulation shows that as the dynamics progresses, the order parameter 
in each horizontal layer converges to nearly constant values essentially deter- 
mined by the temperature assigned to that layer (see Eq. (|55|) ). though some- 
what affected also by elastic effects induced by coupling to n. This is illustrated 
in Figure [5J which shows that the order parameter generally assume smaller 
values at the top and larger values on the bottom, but also varies (slightly) 
within each layer. 

We also study the dynamics of nematic director n. In Figure [31 the nematic 
director on the top surface of the sample is compared at the initial and equilib- 
rium states. For the equilibrium state, three perspectives are given from which 
one can easily discern the evolution of the nematic director. Similar director 
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distributions on the other layers of the sample can be observed. 

The dynamics of the second simulation can be explored similarly. In this 
simulation, all initial conditions and spatial temperature distributions are as in 
the first example, except that one lateral surface of the sample is fixed, and 
the dimensions of the LCE sample are now 4 x 8 x 1 , which is narrower in the 
a;— direction. 

Figure |4] shows the deformation process. Again, the result is very similar to 
that observed in actual experiment (see Fig. 2 of [13]) with the sample bending 
upwards at its free end. As the in first simulation, the deformation proceeds 
through three stages, with the sample also contracting along the initial nematic 
direction and extending in the other two orthogonal directions. This is illus- 
trated the Figs. Ilf; and F. All these phenomena share the underlying physics 
as in the first simulation. 

In Figure [SJ the order parameter distributions within the top, middle, and 
bottom layers of the sample in the equilibrium state are shown. The order 
parameter within each layer is nearly constant within each layer, though with 
small oscillations caused by boundary effects. However, when compared with 
Fig. [2] for the first simulation, one finds that the basic deformation pattern is 
different and asymmetric due to the anchoring boundary condition used in this 
experiment. 

Figure [HI shows the disposition of the nematic director on the top surface of 
the LCEs sample, again comparing the initial and the equilibrium states. For 
the equilibrium state, three perspectives are shown. We see that the nematic 
director bends upwards, and the bending increases with distance from the fixed 
side of the sample. Similar director dynamics are observed in the other layers. 

4 Conclusions 

In this paper, we derived the equations of motion for an LCE sample in the 
long- wave limit, and implemented their numerical solution. Our numerical ex- 
periments demonstrate that the model is capable of describing the dynamics of 
nematic LCEs when exposed to external stimuli such as illumination. 
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Figure 1: The shape evolution of the LCEs sample due to nonhomogeneous changes 
in temperature. Figure A shows the initial state of the LCEs sample, while Figures B, 
C are two intermediate states and Figure D represents the equilibrium state. Figures 
E and F present the shape of the LCEs ^simple at the equilibrium state (Figure D) 
from two different perspectives. In this experiment, the temperature drops linearly 
from the top of the LCEs sample to its bottom while distributes uniformly on each 
horizontal slices, and the temperature spatial distribution is preserved during the 
evolution process. This numerical experiment simulates the real one shown in Figure 
4 in \m. 
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Figure 2: The order parameter (S) distribution for the top (A), middle (B) and 
bottom (C) horizontal slices of the LCEs sample at the equilibrium state. The order 
parameter is close to zero on the top slice while it is close to one on the bottom slice. 
It is the nonhomogeneous distribution of the order parameter that leads to internal 
stress, and thus results in the shape changes of the LCEs sample. Moreover, the order 
parameter slightly varies on each of these slices, suggesting the elastic effect on the 
order parameter. 
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Figure 3: The nematic direction (n) distribution for tlie top slice of the LCEs sample 
at the initial and the equilibrium state. The plot (A) represents the nematic direction 
on the top slice at the initial state. The plots (B), (C) and (D) illustrate the nematic 
direction on the top slice for different perspectives at the equilibrium state. 
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Figure 4: The shape evolution of the LCEs sample due to nonhomogeneous changes 
in temperature. Figure A shows the initial state of the LCEs sample, while Figures B, 
C are two intermediate states and Figure D represents the equilibrium state. Figures 
E and F present the shape of the LCEs 2Anple at the equilibrium state (Figure D) 
from two different perspectives. This experiment shares the same condition as the 
previous simulation in Figure [T] except that a lateral surface of the LCEs sample is 
fixed. This numerical experiment simulates the real one shown in Figure 2 in [13| . 
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Figure 5: The order parameter (S) distribution for the top (A), middle (B) and 
bottom (C) horizontal slices of the LCEs sample at the equilibrium state. The order 
parameter is close to zero on the top slice while it is close to one on the bottom slice. 
It is the nonhomogeneous distribution of the order parameter that leads to internal 
stress, and thus results in the shape changes of the LCEs sample. Moreover, the order 
parameter slightly varies on each of these slices and oscillates near the fixed lateral 
surface of the LCEs sample, suggesting both the elastic effect and the effect due to 
the anchored surface on the order parameter. 
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Figure 6: The nematic direction (n) distribution for the top slice of the LCEs sample 
at the initial and the equilibrium state. The plot (A) represents the nematic direction 
on the top slice at the initial state. The plots (B), (C) and (D) illustrate the nematic 
direction on the top slice for different perspectives at the equilibrium state. 



